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Abstract. The temporal evolution of the entanglement between two qubits evolving 
by random interactions is studied analytically and numerically. Two different types of 
randomness are investigated. Firstly we analyze an ensemble of systems with randomly 
chosen but time-independent interaction Hamiltonians. Secondly we consider the 
case of a temporally fluctuating Hamiltonian, where the unitary evolution can be 
understood as a random walk on the SU{4:) group manifold. As a by-product we 
compute the metric tensor and its inverse as well as the Laplace-Beltrami for S'C/(4). 


1. Introduction 

If two initially separable quantum systems are exposed to random interactions they are 
expected to become entangled, exhibiting random quantum correlations. How do these 
quantum correlations arise as a function of time? To address this question we study the 
entanglement between two qubits subjected to random interactions as a function of time. 
The study of entanglement dynamics under random environments has attracted much 
interest recently, for instance, the emerging entanglement between coupled quantum 
systems through a bosonic heat bath [^. Although our system is oversimplihed in 
comparison with these dissipative systems, we believe that our study may give the 
upper bound for the entanglement under the strong random interactions. 

In what follows we assume that the two-qubit system is initially prepared in a 
well-dehned pure state. As examples we consider two different initial states, namely, a 
non-entangled pure state 


p(0) = |ll){ll| 


( 1 ) 


(a) quenched 


(b) temporal 


Figure 1. Over-simplified cartoon of trajectories on the 5C/(4) group manifold, 
visualizing the two types of randomness discussed in the present work (see 
text). The north pole (red arrow) stands for the identical transformation. 


and in a fnlly entangled Bell state of the form 


p(0) = 10) (01 , 


W = T(|oo) + |ii)), 


( 2 ) 


where {|00), |01), |10), |11)} denotes the canonical qnbit configuration basis. Starting 
with the given initial state the system then evolves unitarily as 


p(i) = c(«)p(0)V(i), 

where the time evolution operator is determined by 17(0) = 1 and 


( 3 ) 


ihdtUit) = Hit) Uit) 


( 4 ) 


with a randomly chosen interaction Hamiltonian. 

Throughout this paper we consider two different types of randomness, namely 

(a) Quenched randomness, where Hit) = H is time-independent. In this case one 
considers an ensemble of two-qubit systems starting from the same initial state, 
where each member evolves by a different but temporally constant Hamiltonian 
drawn from an S'f/(4)-invariant distribution. 

(b) Temporal randomness, where the dynamical evolution is generated by a time- 
dependent Hamiltonian 77(f) which fluctuates randomly as a function of time [^. 
On a single system the resulting temporal evolution of the state vector can be 
understood as a unitary random walk in C^. 

The difference between the two cases is visualized in Fig. In this figure the big red 
sphere stands symbolically for the 15-dimensional group manifold of S'f/(4). Each of 
point on the sphere represents a certain unitary transformation 7/(f) acting on the 4- 
dimensional Hilbert space of the two-qubit system. Starting with t/(0) = 1, which may 
be located e.g. at the ‘north pole’ of the sphere, the temporal evolution t/(f) can be 
represented by a certain trajectory (blue line) on the group manifold. 
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Let us now think of an ensemble of such systems, represented by a set of statistically 
independent trajectories. In the quenched case (a), where a random Hamiltonian is 
chosen at t = 0 and kept constant during the temporal evolution, these trajectories 
are straight, advancing at different pace and pointing in different directions, while in 
case (b) they may be thought of as random walks on the group manifold. At a given final 
observation time the trajectories of the ensemble terminate in different points marked 
by small blue bullets in the figure, each of them representing a unitary transformation. 
Applying this transformation to a pure initial state one obtains a final pure state with 
a certain individual entanglement. In the sequel we are interested in the statistical 
distribution of these final states and their entanglement. 

To quantify the entanglement we use two different entanglement measures. For 
a pure state the entanglement is defined as the von-Neumann entropy of the reduced 
states 


E{t) = -Tr [pi(f)lnpi(t)] = - Tr [psW lnp 2 (t)] , (5) 

where pi, 2 (t) = Tr 2 ,i p(t) denotes the time-dependent reduced density matrix of the 
respective qubit. In cases where the logarithm is too difficult to evaluate we resort to 
the so-called linear entropy 

L(t) = 1 - IV [pl(t)] = 1 - TV [plit)] (6) 

as an alternative entanglement measure. Note that both measures can be obtained from 
the more general Tsallis entanglement entropy 

( 7 ) 

in the limit q ^ and g —)■ 2, respectively. 

Furtheremore, the Renyi entanglement entropy 

H,(t) = (8) 

is of interest. Also this entropy measure generalizes the von-Neuann entropy when 
setting g —1. 

Our main results are the following. In the first case (a), where a temporally constant 
Hamiltonian is randomly chosen, the mean entanglement is expected saturate at a 
certain value for f —)■ cx). Our findings confirm this expectation, but surprisingly we 
observe that the average entanglement overshoots, i.e., it first increases, then reaches 
a local maximum, then slightly decreases again before it finally saturates at some 
stationary value. 

In the second case (b), where the Hamiltonian changes randomly as a function 
of time, the average entanglement saturates as well, although generally at a different 
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level. This saturation level, which is the average entanglement of a unitarily invariant 
distribution of 2-qubit states, has been computed previously in Refs. [5]-|^. Here we 
investigate the actual temporal behavior of the entanglement before it reaches this 
plateau. As a by-product, we compute the metric tensor and its inverse on the 
SU (4) group manifold as well as the corresponding Laplace-Beltrami operator (see 
supplemental material), which to our knowledge have not been published before. 


2. Random unitary transformations in four dimensions 


2.1. Representation of SU{4) transformations 


In what follows we use a particular representation of the group SU{4i) which was 
originally introduced by Tilma et al in Ref. [^. As reviewed in the Appendix, the 
group elements are generated by 15 Gell-Mann matrices Ai,...,Ai 5 , allowing one to 
parametrize unitary transformations U G SU{4:) by 


^ _ giAaai giA2«2 g^Asoa ^iX^a^ giAio«6 g^Asoy ^iX2as 

^ giAaag giAs^io giAaaii giA2«i2 giAaaia gjAgai4 giAisois 


(9) 


where the 15 Euler-like angles a = {oi,..., ois} vary in certain ranges specified in (A.4). 
Applying such a unitary transformation to the non-entangled initial state p(0) = |11) (11| 
one obtains the density matrix 


p(a) = f4p(0)tft 


( 10 ) 


with the components 

Pii(q:) = cos^ ( 02 ) cos^ ( 04 ) sin^ (oe) 

Pi2(ck) = — cos^ (04) sin (202) sin^ (og) 

Pi 3 (q:) = — ^Qg gj,^ ( 204 ) sin^ (og) 

Pi4(q:) = cos (02) cos (04) cos (og) sin (og) 

^ 22 ( 0 :) = cos^ ( 04 ) sin^ ( 0 ( 2 ) sin^ (og) 

^23(0:) = cos (04) sin (02) sin (04) sin^ (og) 

P24{ct) = — e-*(«i-«3-«5) gQg gQg gj,^ gjj^ 

P 33 (q:) = sin^ ( 04 ) sin^ (og) 

^34(0:) = — e*"® cos (cKg) sin (0:4) sin (ag) 

^ 44 ( 0 :) = cos^ (ag) . ( 11 ) 

Remarkably, for p(0) = |11)(11| this density matrix depends only on six angles ai,..., ag 
out of 15. Since the observables investigated in this paper are invariant under local 
unitary transformations, any pure separable initial state will give the same result. 
Therefore, without loss of generality we can restrict ourselves to p(0) = |11)(11|, taking 
advantage of the low number of parameters in this particular case. 
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2.2. Computing averages on the SU(4)-manifold 


In the following section we will consider an ensemble of trajectories of nnitary 
transformations generated by random interactions. Using the above representation, 
each trajectory can be parametrized by a time-dependent vector of Euler angles cx.{t). 
A statistical ensemble of trajectories is therefore characterized by a probability density 
p{cx.,t) to hnd a unitary transformation with the Euler angles ck at a given time t. 

For a given probability density p{cx., t) one can compute the ensemble average of any 
function /(ck) (such as the density matrix p{oi) or the entanglement E{(x)) by integrating 
over the complete parameter space of the S'f/(4) manifold weighted by p{cx.,t)'- 

(fit)) = [ picxC) ficy)dVsu{i)- (12) 

\ hsc/(4) Jvsu(4) 

Here Vsu{ 4 ) is the integrated group volume which serves as a normalization factor while 

15 

dVsu(A) = /i(a) (13) 

i=i 

denotes the volume element on the S'f/(4) group manifold. The actual integration 
measure is defined by the function /x(ck) which depends on the chosen representation. 
Here we use the uniform measure, also known as Haar measure (^, which is by 
itself invariant under unitary transformations||] In the present case of S'17(4) with the 
parametrization defined above the Haar measure is defined by 

pi{oL) = sin(2a2) sin(Q;4) sin®(a 6 ) sin(2a8) sin^(a!io) sin(2ai2) 008 ^( 04 ) cos(a 6 ) cos(aio). (14) 


The total group volume, first computed by Marinov 11 , is then given by 

r9 


ysu(i) — / dVsu{i) — / dcii ■ ■ ■ / dai5/i(Q:) — 


■y/^TT" 


(15) 


In summary, averages over the SU{4:) manifold weighted by the probability density 
p{a.,t) can be carried out by computing the 15-dimensional integral 




dais pia)p{cx,t) /(a), 


(16) 


with the measure (14) and the integration ranges specified in (A.4). The uniform Haar 
measure corresponds to taking piot,t) = 1- 


The transformed state picx) = 11^ p(0) tU is still pure but generally entangled. 
Being interested in the average entanglement of states generated by random unitary 
transformations, it makes a difference whether the entanglement is computed before 
taking average over a or vice versa, as will be discussed in the following. 


f For example, in spherical coordinates the invariant measure of the rotational group SO(3) would be 
given by dlAoji) = g(r, 9)dr d9 df with g{r, (p, 9) = sin 9. 
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2.3. Average of the entanglement 


Let us first discuss the case of computing the entanglement before taking the average 
over all ck. In this case one has to compute the reduced density matrix of the hrst qubit 
for given p{a), defined as the partial trace 

a{cx) = Tr 2 [p(Q:)]. (17) 


For the initial state p(0) = |11)(11| this is a 2 x 2-matrix with the elements 
crii(Q:) = cos^ (04) sin^ (ag) 

<^ 12 ( 0 :) = — sin ( 204 ) sin^ (og) cos ( 02 ) 

- sin ( 02 ) sin ( 2 a 6 ) cos ( 04 ) 

( 722 ( 0 :) = cos^ (cKg) -|- sin^ (q! 4) sin^ (ag) . (18) 

In general the reduced density matrix ^{cx.) is no longer pure and its von-Neumann 
entropy quantifies the entanglement E{cx.) between the two qubits. In order to compute 
the entropy we determine the eigenvalues of a which are given by 


Ki 2 ( 0 :) = = - ± — (256 sin ( 20 ( 2 ) sin (q; 4 ) sin^ (ag) cos^ ( 04 ) cos ( 2 q!i — a^) cos ((Ug) 
’ 2 16 V 

—24 sin^ (ofg) cos ( 20 : 2 ) + cos (2(ag) (8 — 40 sin^ (ofg) cos ( 20 ( 2 )) 
—32 sin^ (2Q;g) cos^ ((U 2 ) cos ( 20 ( 4 ) 

\ 1/2 

-|-32 sin^ (q! 2 ) sin^ (ag) cos ( 4 q! 4 ) -|- 6 cos (4ag) -|- 50 j 
Having determined these eigenvalues, the entanglement of p{cx) is given by 


(19) 


E{a) = — Kj In 


Ki 


1=1 


( 20 ) 


Finally, the entanglement has to be averaged over all trajectories (see Eq. (52)), i.e. 


E = {E{cx))^ 


( 21 ) 


However, if the average of the von-Neumann entropy is too difficult to compute, we will 
also use the linear entropy 

2 

L(a) = £ 2 ( 0 ) = (22) 

i=l 

as an alternative entanglement measure. 


2.4. Entanglement of the average 

Alternatively, we may first compute the average density matrix p = {p{cx))a and 
then determine the entanglement of the resulting mixed state. To this end a suitable 
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entanglement measure is needed. An interesting quantity in this context is Wootters 
concurrence defined by 

C(p) = max(0, Ai - A 2 - A 3 - A 4 ), (23) 

where Ai,..., A 4 are the decreasingly sorted square roots of the eigenvalues of the matrix 

A = p(cr^ ® ® (J^) . (24) 

In this expression is the Pauli matrix while p* denotes the complex conjugate of a 
without taking the transpose. 

From the concurrence one can easily compute the entanglement of formation of the 
mixed state, which is given by 

E^(p) =-61n6-(l-6)ln(l-6), (25) 

where b=\ + \- C{pY. 


3. Quenched random interactions 


In the case (a) of quenched randomness each element of the ensemble is associated with 
a time-independent random Hamiltonian H. Since the spectral decomposition 

4 

= (26) 

i=i 

of a randomly chosen Hamiltonian is always non-degenerate, the time evolution operator 
can be written as 

4 

= (27) 

i=i 

Hence the state of the system evolves as 


Pit) = C/(«)p(0)C/t(i) = (28) 

j,k=l 

where p( 0 ) denotes the initial state. 


The Hamiltonian itself has to be drawn from a certain probabilistic ensemble of 
Hermitian random matrices [2,13 . Here the most natural choice is again the Gaussian 
unitary ensemble (GUE). This ensemble has the nice property that the probability 
distributions for the eigenvalues Ej and the eigenvectors | 0 j) factorize and thus can be 
treated independently. More specifically, the eigenvalues are known to be distributed as 


P(Ei,...,E4) oc 

n>m 


(29) 
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where A = ^ is a constant determining the width of the energy flnctnations and 
therewith the time scale of the temporal evolntion. In the following the corresponding 
average over the energies will be denoted by On the other hand, the orthonormal 

set of eigenvectors is randomly oriented in the fonr-dimensional Hilbert space according 
to Haar measnre, independent of the eigenvalnes. If one defines the qnbit basis 

{|1>,|2),|3>,|4>}:={|00),|01),|10>,|11>} (30) 

this average can be carried ont by setting 

I'^j) := Ul\j) (31) 

and integrating over all Enler angles ct according to the Haar measnre (see Appendix B). 
This average will be denoted by . .)a- The total GUE average thns factorizes as 

(...)gue = (32) 


3.1. Entanglement of the averaged density matrix 


Let ns now compnte the average density matrix 


S 




/ GUE \ / 

j,k=l ^ - 

^jk 


E 




(33) 


T 


jk 


First we compnte the average over the energies 

1 ^+oo 

Rjk = J^j dEi---dE4PGUE(Ei,...,E4)e-*('^^-'^'')*, (34) 

where Af = dEi ■ ■ ■ dE 4 Pgue(-Ei, • • •, P 4 ) = ^ is the normalization factor. This 


2 A8 

1 for j = k 

/(r) for j y^k 


leads ns to the result 

Rjk = f{r)+ (1 -/(r)) 5 jfc = 

where we defined the scaled time 
T := t/\/2A 
and the function 

/(r) = (-2r^° + 25r® - 128x6 + 276r^ - 288x2 + 72) . 


(35) 


(36) 


(37) 


Thus Eq. (33) reduces to 

[pit)) 


/ GUE 


4 4 

/WEtj<.+ (i-/W)Et« 

j,k=l j=l 


( 38 ) 






What remains is to determine the operators 


= {\4>,){<h\m\4>t){4>k\) = (uimwo. mu}.\kmuo) . ( 39 ) 


Obviously, the first sum in Eq. (38) is given by 


4 

= (t 2 u.p( 0 )irlu,) = (p( 0 )) = p( 0 ). 


(40) 


j,fc=l 


As for the second sum in Eq. (38), we note that the distribution of eigenvectors 
\4’i) = U}^\i) is invariant under a permutation of the basis vectors \j), hence the 

four operators Tjj coincide. Moreover, under a unitary transformation V G S'f/(4) they 
transform as 


= (vui\]){]\Uc p( 0 ) 


= (ul\j)(j\Uc (v"p(0)V) l2|j)(j|£/< 


(41) 


where we have used the invariance of the GUE-eigenvectors under the replacement 
Ua —t UaV. This means that Tjj is invariant under V if and only if V commutes with 
the initial state p(0). For a pure initial state this implies that Tjj has to be a linear 
combination of the identity and the initial state itself, i.e. Tjj = al + bp{0). The linear 
coefficients a and b can be determined as follows. On the one hand, the identity 


1 = Tr 


(P(^))c 


( | 38 | 40 | | 


/(T) + (l-/(T))Tr[ET 


rj 


(42) 


implies that Tr ~ hence 4a + 6 = 1/4. On the other hand, we note that 


Tr 


p(0)T 


33 


= (Tr 


P(O)|0i)(0j|p(O)|0j)(0j 


= (0.|p(O)|0,y (43) 


is invariant under unitary transformations of p(0) and independent of j, hence we may 
choose j = 4 and p(0) = |4)(4| to obtain 

lv[p(0)T,,] = = (cos‘‘(a6))_^ = i, (44) 

giving a = b = Therefore, we arrive at the convex combination of 1/4 and p(0) 


m-{m) = izi(i),^i±iE) 


p(0) 


(45) 


with /(r) given in Eq. (37) and t = t/ \/^, which holds for any pure initial state p(0). 
As expected, the averaged state lies on the segment between the initial state p(0) and 
the maximally mixed state 1/4 due to the symmetries of the Haar measure of SU (4) . 

Having computed the mixed state of the ensemble p{t) we can now compute the 
corresponding entanglement of formation as a function of time. For a non-entangled 
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— EF(0)=ln(2) 

- E/=(0)=0.520 

- Ef(0)=0.344 

- Ef(0)=0.173 

— Ep(0)=0 


Figure 2. Quenched case: Analytically calculated entanglement of formation 
of the averaged density matrix in Eq. (45) as a function of the scaled time 
r = t/\/2A using different initial conditions p(0) = |V’)(V'| of pure states 
I'0c) = (c, 0, 0, a/I — c^) with c = {0,0.204, 0.330, 0.464, l/\/2} from bottom 
to top. In particular, the purple and the blue line represent a non-entangled 
and maximally entangled initial state, respectively. 


initial pure state p(0) = |11)(11| we find that Ep^pit)) = 0 for all times. However, if we 
start from the Bell state (|^ with the initial entanglement Ep{p{fi)) = ln2 0.693 we 
find numerically that the entanglement first decreases and vanishes at the finite scaled 
time T 0.4997. (see Fig. |^. 


3.2. Average of the individual entanglement 


Instead of computing the entanglement of the average density matrix, let us now 
compute the average of the individual entanglement of each trajectory, i.e. the 
entanglement is computed before taking the GUE average. Although the von-Neumann 
entanglement entropy of the individual pure states would be straight-forward to compute 


(see (20)), we did not succeed to compute the average. For this reason let us consider 
the GUE average of the linear entropy 


(m) =i-{pi{t)) 

\ / GUE \ / ( 


(46) 


where pi{t) denotes the reduced density matrix of the first qubit. In the qubit basis (30) 
this can be rewritten as 


Lit) 


= 1 


( (E/^ldWl7/5) 


/4,/3,7d=l 


(47) 
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Inserting (28) and exploiting again that the GUE average factorizes, we get 


( 48 ) 


'm / ’ 


(cf* {(l)j\p{0Mk) cf cf* (7|p(O)|0™)c^ 

\ I UL 

where = {(f)j\fi(3). For the initially non-entangled state p(0) = |11)(11| this 
expression reduces further to 

2 4 

^L(t)'^ = 1 — ^ ^ ^ ^^-iiEj-Ek+Ei-Em}TV^'^ 

//,/ 3 , 7 ,( 5=1 j,k,l,7n=l 

Rijkl ) 


(49) 


^ ^ S S ‘-'I ^rn j ^ 


jklm 


jklm 

with the scaled time r = t/y/^A. As shown in Appendix D, the averages Rijki{T) 
can be computed directly by integration over the given probability distributions in GUE, 
leading us to the final result 

(L(r))onE = (32rS - 128r'^ + 168r" - 72^ + 9) (50) 

1 


(-2r^° + 25r® - 128r® + 276r^ - 288t^ + 72) 
(-54r^° + 387r® - 832r® + 828r^ - 288r^ + 24) 


420 
1 


1 s 

(-256r^° + 800r® - 1024r® + 552r^ - 144r^ + 9) + — 


315 


70 


This function is plotted in Fig. As one can see, the linear entanglement entropy 
(black line) first increases rapidly, then reaches a local maximum (L) ~ 0.199936 < 0.2 
at Tmax — 0.817377, then decreases again and finally saturates at the value 


lim (L(r))GUE = 13/70 ~ 0.1857. (51) 

T^OO 

Because it would need much more effort to calculate the analytical the linear entropy 
for different initial states analytically, we used numerical methods. The results are 
compared in Fig. As one can see clearly, all the lines tend to touch the limit value at 
the fixed time r ~ 0.817 and the curves do not intersect. 

Fig. [^explains the meaning of this result: Each single trajectory of the ensemble 
on the surface is deterministic with a given initial point, direction and velocity. Since all 
members of the ensemble share the same initial starting point on the upper half of the 
sphere, the probability for finding the walkers can be slightly higher on the upper half in 
the long-time limit because all trajectories will periodically return to this point. This is 
why the limit depends on the initial state and therefore deviates from the Haar measure. 
The bump can be seen as a transient state in which the probability distribution seems 
to be almost randomly distributed before saturating. 
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. L(7-)=1/5 

-j-=0.817 

— Analytic 

— Numeric 1 

— Numeric 2 

— Numeric 3 


r 


Figure 3. Average linear entanglement entropy according to Eq. ( |50[ ) as a 
function of the scaled time r = tl\f2A with c = {0,0.204,0.464, l/\/2} from 
bottom to top. The numerical data (see. Appendix E) has computation 


errors smaller than the thickness of the lines. The upper right panel shows 
the magnified area around the touching point at (marked by the dashed 
line). 


Note that in contrast to the case discussed before (see Fig.[^ the system remembers 
its initial state, saturating at different levels of entanglement in the limit t —)■ oo. 


4. Time-dependent random interactions 


Let us now consider the case (b) of a temporally varying Hamiltonian, where the state 
vectors of the ensemble perform a unitary random walk on the S'[/(4) manifold. In 
this case the quantity of interest is the probability distribution p{a., t) to hnd the time 
evolution operator U{t) with the Euler angles cx at time t. This probability distribution 
allows one to compute the ensemble average of any function /(ck) (such as the density 
matrix p{(y.) or the entanglement Ect) by integration over the complete S'17(4) volume 
weighted by p{cx,t), i.e. we have to compute the integral 


fit)) = 


hs{/(4) Jvsu(4) 


f{cx)p{cx,t) dVsu{A) 


(52) 


over the ranges specihed in (A.4), where (iVsu{A) denotes the volume element according 
to the Haar measure dehned in Eq. (13). 
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4-1. Expected average entanglement of a uniform distribution 

Before studying the temporal evolution in detail, let us consider the limit t ^ oo, 
where we expect the state vectors to be uniformly distributed on the group manifold. 
Since such an ensemble is by itself invariant under unitary transformations, the state 
vectors are distributed according to a Gaussian Unitary Ensemble (GUE). Starting 
from this observation. Page conjectured a closed expression for the expected average 
entanglement of an arbitrary bipartite quantum system with Hilbert space dimensions m 
and n, which was later proven rigorously by Foong, Kanno, and Sen . This formula 
describes the average entanglement of a random pure state between two subsystems 
with Hilbert space dimensions m and n: 



Applying this formula to a random two-qubit system, one obtains 

= {E,,,) = i. (54) 

This is the average entanglement between two qubits in a randomly chosen pure state. 


4-2. Heat conduction eguation on the S'f/(4) manifold 


In order to compute the probability distribution p{a,t) one has to solve the heat 
conduction equation 


dpicx,t) 

dt 


DAaP{ct, t) = 0 


(55) 


on the curved SU (4) group manifold. Here D denotes the diffusion constant while Aq, is 
the so-called Laplace-Beltrami operator which generalizes the ordinary Laplacian on a 
curved space. As stated by this Laplace-Beltrami operator is the Markov generator 
of the unitary brownian motion (see Appendix F). On a Riemannian manifold with the 


metric tensor gij the Laplace-Beltrami-Operator A is given by 

1 


A/ = 


yisi 


d, V\a\9%f 


(56) 


where g^^ are the components of the inverse metric tensor and \/\^\ = \/det g. 

To our best knowledge the explicit expressions for gij, and Aq, have not been 
published before. This is perhaps due to the fact that the formulas are so complex that 
even powerful computer algebra systems such as Mathematica® are not able to compute 
the inverse metric directly. Instead one has to invert the matrix manually element by 
element. Our explicit results are included in the supplemental material attached to this 
paper. 
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4-3. Early-time expansion 


The solution p{oL^t) of the heat conduction equation (55) and as well the averaged 
function {fit)) can be expanded as a Taylor series around t = 0: 


lY) ‘~\TI 

p{a,t) = 


71=0 


(67) 


t=0 




71=0 


(58) 


i=0 


Using (52) we can compare the coefficients of the both Taylor series, giving 
/ \ 1 


dt"^' 


t=o 


dVsu{4) /(a) ^p{a,t) 


t=o 


(59) 


where dVsu{ 4 ) = V^Md^^a denotes the volume element defined in Eq. (13). Using the 
heat equation 


we can replace the partial derivative, obtaining 

/ „, ^\ _ 1 

t=o Vsu{4) 


&^ / \ 
wMf 


D- / dUsn(4) /(a)A>(a,f) 


t=o 


(60) 


The r.h.s. is an integral over derivatives of the probability density p{a,t) evaluated at 
t = 0. If all trajectories start at ckq it is easy to see that this probability density at 
t = 0 is given by 


p{(x,t = 0 ) = - (Xq). 

V\9\ 


(61) 


Inserting this expression into ([60j) the integral can be evaluated by partial integration, 
giving 


dn 

dt'^ 


{fit)) 


t=o 


D-A-f{cx) 


Q:=Q :0 


(62) 


f-f- Average of the density matrix 


Using (10) we calculate the derivatives A2 p(q:). We find that 

A^p(a) =-8p(a)+ 2-1. (63) 

Therefore we can express all higher derivatives of n > 1 in terms of the hrst derivativ^ 
(9" / \ 




= D^{-8)^-^AMa) . 


t=o 


§ Although this remarkable property calls for a deeper reason, we have no convincing explanation so 
far. 
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Hence, the solution of the averaged density matrix can be written as 




( 65 ) 


Using the non-entangled initial state |'0o) = |11) (which corresponds to taking ckq = 0) 
this result specializes to 


(pit)) 

\ / Q;o=' 



V 


l„-8Dt 

0 

0 

0 


1 

4 


0 

1 „-8Dt 

0 

0 


1 

4 


0 

0 

Ip-SDt 

A^ 

0 


0 \ 

0 

0 


1 I 3 -8Dt 

4 t" 4^ 


.( 66 ) 


As can be seen, this density matrix relaxes exponentially and becomes fully mixed in 
the limite t —>■ oo. 


4-5. Average of the linear entropy 


The same calculation can be carried out for the linear entropy defined in (22). Here we 
find that 


AiL(a) =-20L(a)+4. (67) 

For this reason, the calculation is completely analogous, giving 

{m) = t + (^L(q„) - (68) 

The result for an non-entangled initial state L{cx.q) = 0 is therefore 

(m) = t (69) 

For a fully entangled initial state L{cx.o) = 1/2 we get 

(m) = ^ + ( 70 ) 


In both cases the averaged linear entropy tends towards | which coincides with the value 
of a randomly distributed density matrix according to Haar-measure (cf. Ref. (5|[l4]). 

4 . 6 . Average of the Tsallis entropy 

Because of computational difficulties we calculated the Tasllis and von-Neumann entropy 
only calculated to first order in t. Moreover, we restricted ourselves to the case of a 
fully entangled initial state since the Taylor expansion fails for a non-entangled initial 
state due to a logarithmic factor in time at t = 0 . 
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T 


. Limit: E(t)=1/3 

— E(0)=0, unentangled 

E(0)=log(2), fully entangled 
E(0)=0.1731 

- E(0)=0.5204 

- E(0)=0.2996 

- E(0)=0.3333 


Figure 4. Numerical results (see Appendix E) of the von-Neumann entropy 
compared to the analytical approximations (dashed), starting from pure initial 
states p(0) = \ipc){ipc\ with c = {0, 0.203, 0.300,0.322,0.464, l/\/2} (see Fig. [^. 
The error of the numerical calculation is smaller than the thickness of the lines. 
The scaled time is defined as r = 2Dt. 


Using the definition of the Tsallis entropy ([^ and following the same lines as 
outlined above, one ends up with the first order approximation 


E,{t)) = 


1 - 2 ^-'' 




(71) 


In the limit g —)■ 2 this reproduces the hrst order terms of (70) whereas for g —>■ 1 we 

(72) 


obtain the hrst-order approximation of the von-Neumann entropy 

=\og2-QDt + 0(f). 


By direct approximation on can derive the hrst order in r = 2Dt of the von-Neumann 
entropy directly for a few timesteps using an initial state l-^o) = (cos <)), 0 , 0 , sin(/)) with 
0 e (0,7r/4): 

(^(^)) = (- log(sin(0)) - log(cos(0)) - cos(20) log(cot(0))) -h (73) 

(2 cos(4(/)) sec(20) log(cot(0)) — 1) ■ r -|- 0{t^) 

For a fully entangled initial state (0 —>■ 7 r/ 4 ) this reproduces Eq. ( |7^ . Unfortunately, 
the expansion does not converge for 0 —)■ 0 , which is why one has approximate an 
unentangled initial state directly, obtaining 

<^E(r)) = r (7 - logr)-k C>(r^), (74) 

where 7 is Eulers constant with numerical value of 7 ~ 0.577. 

To clarify these results Fig. shows numerical calulations and the analytical hrst- 
order approximations as dashed lines for several initial states. 
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It is interesting to note that the von-Neumann entropy shows a bump if one starts 
near the limit value of 1/3, whereas the linear entropy does not have such a behaviour. 
But since one can interpret the linear entropy as lowest-order expansion term of the 
von-Neumann entropy around a pure state, it is not surprising that the von-Neumann 
entropy shows a more complex behavior. Moreover, a fixed initial state with a von- 
Neumann entropy of 1/3 is obviously not yet randomly distributed in 517(4). 


J^.l. Average of the Renyi entropy 


We perform similar calculations as above for the Renyi entropy as defined in (|^. This 
leads to the following approximation in first order for small times using a fully entangled 
initial state: 


(75) 


Hq{t) ) = log2 — 6Dqt + 0{r 


For g —)■ 1 we obtain the von-Neumann entropy as already computed in Eq. (72) 


Setting q = 2 one can compute the first order approxiation first order in r = 2Dt 
of the Renyi entropy for the state l-^o) = (cos(/), 0, 0,sin(/)) as defined above: 

<^772(r)^ = log4 - log(cos(40) -F 3) (76) 

(28 cos(40) + 3 cos(80) + 1) _ 25 

(cos(40) -|- 3)2 

Note that this equation holds for all (j) G (0,7r/4). For a fully entangled initial state 


(0 —>■ 7r/4) this reproduces Eq. (75) using q = 2, and for an unentangled initial state 


(0 —)■ 0) one receives 

(// 2 (r)) =2r+0(T") 


(77) 


5. Discussion 


In this paper we have studied how the entanglement of a two-qubit system subjected 
to random interactions evolves in time. We have considered two different types of 
randomness, namely, quenched and time-dependent random interactions, starting either 
from a separable or from a maximally entangled initial state. The main results are the 
following: 

• Quenched random interactions: Since entanglement measures are non-linear, 
it makes a difference whether the average over the quenched disorder is carried out 
before or after evaluating the entanglement |j|] In the first case we can compute the 

II In ffeneral, it would be interesting to see whether there exists an inequality relating the quantities 
{Eip{t)))andE{{p{t))). 
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averaged density matrix explicitly (see Eq. (45)), finding that the entanglement 
of formation decreases monotononsly, see Figj^ In the second case we can only 
compnte the so-called linear entropy (50), which is fonnd to overshoot before it 
satnrates. 


• Time-dependent random interactions: As ontlined in the Introdnction, this 
problem is eqnivalent to solving a random walk on the 517(4) gronp manifold. We 
find exact expressions for both the averaged density matrix and the averaged linear 
entropy, confirming that these qnantities vary exponentially with time. For the 
averaged von-Nenmann entropy as a special case of Tsallis entropy, however, we 
conld only obtain a first-order approximation. 


It shonld be noted that taking the average after evalnating the entanglement leads 
to resnlts which are not directly accessible in experiments. The reason is that the 
compntation of entanglement in pnre-state systems reqnires the knowledge of the fnll 
density matrix in one of the snbsystems. This matrix can only be measnred by means of 
repeated experiments nnder identical conditions, which in onr case wonld mean to nse 
the same realization of randomness. Having estimated the entanglement this resnlt has 
then to be averaged over different realizations of randomness. In experiments, where 
the randomness differs npon repetition, one wonld instead obtain a fnlly mixed density 
matrix withont any information abont entanglement. 


As a by-prodnct, when solving the random walk problem on the 517(4) manifold we 
had to compnte the corresponding Laplace-Beltrami operator which in tnrn reqnired to 
compnte the metric tensor and its inverse. The compntation of the metric was extremely 
difficnlt and to onr best knowledge has not been done before. The resnlting expressions 
are lengthy (see Snpplemental Material). Snrprisingly, we finally arrive at very simple 
resnlts (see Sect. 4.4). This indicates there might be a deeper mathematical strnctnre 
behind the problem that we failed to nnderstand so that the brnte-force calcnlation 
presented here is perhaps not really necessary. It wonld be interesting to investigate 
this point in more details. 


Appendix A. Representation of 517(4) in terms of Euler angles 

The Lie algebra of the symmetry gronp SU (4) is defined in terms of 15 Gell-Mann-like 
generators Ai,..., A 15 which can be represented in varions ways. In this paper we nse 
a representation snggested in j9, 15 . Denoting by Eij a 4 x 4 nnit matrix in which the 
element at position i,j is 1 and all others are zero, this representation is given by 


■1 = 

Ei 2 -|- E21 

A2 = '1(^21 — E12) 

-^3 = 

Ell — E22 

4 = 

Fi 3 -|- E31 

As = i{E 3 i — Elf) 

Ae = 

E23 + E32 

7 = 

i{E‘i2 — E23) 

As = + -^22 — 2E33) 

Ag = 

Eli + Eil 

ao = 

= i{E 4 i — E14) 

All = F24 -|- E42 

Ai2 = 

= i{Ei 2 — E24) 

13 = 

= E43 -|- E34 

Ai 4 = i(E43 — E34) 

-^15 = 

= + E22 


+ -£^33 ~ 3E44) 
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The generators are Hermitian and obey the relations 


Tr[Ai] = 0 , Tr[A-] = 2 , [A^, A^] = ^ fjkiXi , 


i=i 


where fjki are the S'17(4) structure constants defined by 


(A.l) 


fjki ~ [Aj,Afc]Ai 

Group elements U G S'f/(4) can be generated by 

U[ct) = g*^ 2«2 giAsas ^1X3015 giAioae g*A 3«7 ^1X303 

^ ^iXsag giAsaio ^iXgan ^iXgaig giAsaia giAgai4 giAisais 


(A. 2 ) 


(A.3) 


where a = {ai, 02 ,..., 015 } is a set of 15 parameters analogous to Euler angles. Their 
ranges are 


a2, 0:4, 0:6, Oig, aio, 0:12 

(y.\^ 0^7, 0^11 

as, as, ag, ais 

ai4 

ttl5 


e [0,7r/2] 

G [0, tt] 

G [0,27r] 
e [0,^7r] 

G [0,yp7r]. 


(A.4) 


According to [9], the main advantage of this parametrization in these ranges is that it 
provides a coverage of the group manifold without overlaps, i.e., the parametrization 
does not overcount the manifold. There exists also a more recent representation with 
similar properties which is more symmetric and transparent in the definition of the 
angles 16 . However, the complexity of the formulas turns out to be comparable. 


Appendix B. SU(4) Haar measure 

A measure /i on a compact group G is called Haar measure if it is translation-invariant 
under the group itself, i.e., for any subset S' C G we have ^{g o S) = g{S) for all g E G. 
Loosely speaking the Haar measure may be thought of as a constant probability density 
on the group manifold. In a given parametrization this requires to define a suitable 
invariant volume element on the group manifold. In the present case of SU (4) with the 
parametrization defined above this volume element is given by 

dGsc/(4) = /i(Q:) dai dag ■■■ dai5 (B.l) 

with 

= sin ( 2 a 2 ) sin (a 4 ) sin® (ae) sin (2a8) sin® (aio) sin ( 2 ai 2 ) (B.2) 

X cos® (a 4 ) cos (ag) cos (aio) • 
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The total group volume of 517(4), first computed by Marinov [TT], is then given by 

r9 




SU{4) 


= / dVsui4) = / dai-- - / dai5/i(Q:) = 


V27r^ 


(B.3) 


where the integration is carried out over the ranges specified in (A.4). 


Appendix C. Derivation of the metric tensor, its inverse and 
Laplace-Beltrami-Operator of 5f/(4) 


First we derive the metric g of SU{4) using the representation of the manifold given 
by (g. Using the induced scalar product of matrices one can express the inhnitesimal 
line element as 


= Tr [dUdU^] 


(C.l) 


In general the line element on a Riemann manifold is defined as 
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ds^ = Qijdaidaj 
i,j=^ 


(C.2) 


Thus by equating coefficients in (C.l) and (C.2) one obtains the matrix elements of the 


metric g. This metric with coefficients gij has to be inverted to g'^^ in order to calculate 


the Laplace-Beltrami-Operator in (56). To this end one needs the square root of the 
determinant of the metric \/|^, but this is already given by (B.2). The inversion is 
difficult and was done manually for individual matrix elements. 


Appendix D. Integration of Rijkiij) and Tijuxajie 


The following table lists the appearing results of the integration of Rijkiir) for the 
different cases of the indices. 


Averaging Tijkixa^e for the initial state |V’o) = |4) means to average the 4-component of 
the unitary matrix defined by (31). Thus one can rewrite Tijkixajie as 


TijklXa^e — {Ul^Uj^iUl^JJi^JJi^2(X-l)+o-Ul2(^ii_i)j^JJi,2(l3-i)+eUl2(^X-l)+e)oc (D.l) 


where we use the mapping for the qubit basis (30) 


ot) 0 \(d) — |2a -)- /3 -)- 1). 


(D.2) 


Not all of these integrals have to be carried out. Since the unitary matrix is randomly 
distributed according to Haar-measure, the rows and columns can be swapped as desired 
without changing the result, e.g. T 12141222 = 7'i24ii222- 
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Index propertjes 

Rjklm 

^ k ^ 1 ^ m 

(9 - 36r2 + 42r^ - 16r® + 2r^) 

< < < 

(1152 - 2304r2 + 1104r^ - + 25r® - 

{j = 1 A k ^ m) W 
{k = m Aj ^ 1) 

^e-^(384 - 2304r2 + 3312r^ - 1664x6 + 387r® - 27ri°) 

j = 1 A k = m 

ie-2^'(9 - 72x2 + 138x^ - 128x6 + 50x® - 8x^0) 

{j = k Al = m)\J 
{j = m A k = 1) 

1 


Table Dl. Results for the integration of Rijki{T) 


Appendix E. Numeric calculations 


Depending on the case of the randomness we use different approaches for the numerical 
calculation. 


In quenched randomness case it is possible to compute p{t) directly for a single 
member of the ensemble using Eq. (28), because H is random but constant. 


• In the temporal case we have to choose a new random Hamiltonian at each numerical 
time step dt. The unitary operator is chosen as U = A~^ ■ with A = 1 + iH^ 
to ensure a unitary transformation. To ensure scale invariance in time, one has to 
scale a chosen H by l/\/di as expected for noise terms. 


Afterwards one can carry out the average over the ensemble to compute the averaged 
density matrix and the entanglement of formation, or to obtain the averaged 
entanglement measure (linear, von-Neumann). 

The ensemble of the random Hamiltonians H is the GUE. To obtain a H & GUE 
one has to chose a Hermitian matrix containing 

• random numbers z with z G Af{0, on the diagonal entries and 

• complex random numbers with Zi,Z 2 G A/'(0, cr^) on the off-diagonal entries, 

where N'{0,a‘^) is the normal distribution. 


Appendix F. Motivation for the diffusion equation on S'f/(4) 


In the following we would like to sketch how to derive the diffusion equation (55) 
starting from the time evolution operator U{t) in (|^ in the case of temporal random 
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interactions. The time evolution operator U{t) can be written as product of iV = T 
unitary transformations, that act with a new random Hamiltonian Hi in the time step 
i with constant infinitesimal time interval dt: 


IV’W) 


t/(f)|V’(0)) = UN{dt)U^_,{dt)---U,{dt)Uo{dt)\'iP{0)) 




g-iHt/dt 


\t'=0 


\m) = 
I IV'(O)) 



IV'(O)) = 


(F.l) 


with Uo{dt) = 1, so Hq = 0. In this expression we rename the position index i by a 
temporal index f. Moreover, we define 


dHt’ = Hfdt, 


(F.2) 


SO that 

t 

U{t) = (F.3) 

t'=0 

In the following we want to show that the defined time evolution operator U{t) equates 
a stochastic differential equation of a unitary Brownian motion (UBM), it the randomly 
chosen Hamilton Ht' operators are drown out of 

Hf = (F.4) 


Here the — .^gue are randomly chosen matrices of the GUE using the normal distribution 
Ar(0,(T^). The sign can be chosen freely because of the symmetry around zero. Using 
Eq. (|F.2l) and (IF.41) we get 


dHf/ — .^GTiF.V^df . 


(F.5) 


These are random matrices according to the GUE whose components are random 
numbers out of A/'(0, a^df). Introducing a normalized time t ^ the random 
distribution changes to A/'(0, dr), which is why dH^-i generates a Wiener process. 


Having a look at the differential dU(r) we conclude that 


dU (r) = U{t + dr) — U (r) 

( T+dr \ / '^ \ 

n uv) - n u{r') 

t '=0 / \ r '=0 / 

= (^nr(r')j(C/(r + dT)-l) 

= f/(r) - 1) 

= iU{T)dHr - ^f/(r)dr. 


(F.6) 
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where we expanded the exponential fnnction in the last step for dr —>■ 0 using (F.5). 
This is a stochastic differential equation for the Ito process {f/(r)}T-gK+ and therefore 
describes a unitary Brownian motion on S't/(4) (see also 18,1^). 


As stated by 17 the Markov generator A of this diffusion process is given by the 
Laplace-Beltrami-Operator A on the Riemannian manifold A = |A. As described 
above, the metric is induced by the matrix scalar product. Since the Laplace-Beltrami- 
Operator is independent of the choice of the basis, we use an operator Aq; that is 
parametrized by the Euler angles ct, getting 

dp{cx, t) 


dr 


= 2^ap(a,T) 


(F,7) 


Using the rescaled time t = ^ and the definition D = ^ we end up with the diffusion 
equation (55). 
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